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We derive expressions, in terms of "polar shapelets", for the image distortion operations 
associated with weak gravitational lensing. Shear causes galaxy shapes to become elongated, 
and is sensitive to the second derivative of the projected gravitational potential along their line 
of sight; flexion bends galaxy shapes into arcs, and is sensitive to the third derivative. Polar 
shapelets provide a natural representation, in which both shear and flexion transformations are 
compact. Through this tool, we understand progress in several weak lensing methods. We then 
exploit various symmetries of shapelets to construct a range of shear estimators with useful 
properties. Through an analogous investigation, we also explore several flexion estimators. 
In particular, some of the estimators can be measured simultaneously and independently for 
every galaxy, and will provide unique checks for systematics in future weak lensing analyses. 
Using simulated images from the Shear TEsting Programme (STEP), we show that we can 
recover input shears with no significant bias. A complete software package to parametrize 
astronomical images in terms of polar shapelets, and to perform a full weak lensing analysis, 
is available on the world wide web. 
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1 INTRODUCTION 

Weak gravitational lensing is a powerful method to map the dis- 
tribution of mass in the Universe, regardless of its nature or state 
(for reviews see Mellier 1999; Bartelmann & Schneider 2001; Re- 
fregier 2003). The apparent shapes of background galaxies become 
distorted as their light travels near mass concentrations along their 
line of sight to the Earth. The well-known shearing of galaxies, in 
which intrinsically circular sources would be seen as elongated el- 
lipses, is induced by an amount proportional to the second deriva- 
tive of the projected foreground gravitational potential. Such dis- 
tortion has been measured around individual galaxy clusters (e.g. 
Wittman et al. 2001; Wittman et al. 2003; Bacon & Taylor 2003; 
Bradac et al. 2005 ; Wittman et al. 2006) and, in a statistical fashion, 
by large scale structure (recent measurements include Massey et al. 
2004; Van Waerbeke et al. 2005; Heymans et al. 2005; Jarvis et al. 
2006; Hoekstra et al. 2006; Semboloni et al, 2006; Hetterscheidt et 
al. 2006; Schrabback et al. 2006; Kitching et al 2007; Massey et 
al. 2007b). 

A higher-order effect, known as "flexion", is also emerging as 
a probe of the distribution of mass on small scales, and particularly 
in the inner cores of galaxy clusters (Goldberg & Natarajan 2002; 
Irwin & Schmakova 2003; Goldberg & Bacon 2005; Bacon et al. 
2006; Irwin & Schmakova 2006; Okura, Umetsu & Futamase 2006; 
Goldberg & Leonard 2007). Variation in the shear signal across the 



width of a background galaxy causes bending in its apparent shape. 
This is the next term in a lensing expansion that leads towards the 
formation of an arc, as in strong lensing. The flexion is sensitive to 
the third derivative of the projected gravitational potential. 

Precise image analysis techniques are required to detect weak 
gravitational lensing, because the shapes of galaxies are changed 
by the effect by only a few percent. In fact, the lensing contribution 
to the shape is about an order of magnitude smaller than the dis- 
persion of galaxies' intrinsic morphologies and the spurious distor- 
tions introduced by typical imperfections in telescopes. The widely 
used shear measurement method by Kaiser, Squires & Broadhurst 
(KSB, 1995) has been successful in many contexts, but contains 
several documented shortcomings: it is found to be insufficiently 
accurate to measure shears with a desired accuracy of less than 1 % 
(c.f. Bacon et al. 2001; Erben et al. 2001; Heymans et al. 2005 
(STEP1); Massey et al. 2007a (STEP2)), and it is mathematically 
ill-defined for realistic Point Spread Functions (c.f. Kaiser 2000; 
Kuijken 1999; Hirata & Seljak 2003). 

Several new shear measurement methods are being developed, 
to fully exploit future space-based weak lensing surveys with HST 
or the proposed SNAP, DUNE or IDEM missions, and ground- 
based wide-field surveys such as those with Megacam, CTIO DES, 
VISTA darkCAM, Pan-STARRS and LSST. A review of the var- 
ious shear measurement methods is found in STEP2, along with 
their division into "active" and "passive" categories. Active tech- 
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niques work by modelling galaxies as intrinsically circular, then 
shearing the models until they most closely match the observed 
ellipticities. Passive methods work by measuring the apparent el- 
lipticities of objects as well as higher order shape statistics, which 
are used to calibrate the ellipticities. 

Flexion measurement methods are still in relative infancy. 
Initial attempts to mathematically describe the flexion distortion 
(Goldberg & Natarajan 2002; Irwin & Shmakova 2003) were 
formidably complicated. A passive estimator has been constructed 
by Okura, Umetsu & Futamase (2006), and further expanded by 
Goldberg & Leonard (2007). A completely different, probabilis- 
tic approach is taken by Irwin & Shmakova (2006) and Irwin, 
Shmakova & Anderson (2006). However, several important fea- 
tures in these approaches remain to be developed, and they remain 
mathematically complex; it is therefore desirable to find a formal- 
ism which allows maximum physical insight into the problem. An 
advance towards this was made by Goldberg & Bacon (2005), who 
related flexion to the formalism of Cartesian shapelets (Refregier 
2002; Refregier & Bacon 2002). Shapelets contain all the mechan- 
ics necessary to deconvolve galaxies and flexion estimators from 
the effects of a PSF The active method of Goldberg & Bacon 
(2005) and Goldberg & Leonard (2007) has been used to success- 
fully detect the flexion signal. The mathematics has a simpler form, 
although it is still not as elegant as possible. 

Here, we present the image manipulations of lensing theory in 
terms of the "polar shapelets" formalism (Refregier 2003; Massey 
& Refregier 2005). This suggests a complete, orthonormal set of 
basis functions into which any galaxy shapes can be decomposed. 
It also provides a neat way to deconvolve arbitrary galaxy shapes 
from an arbitrarily complicated PSF, so we can set out under the 
assumption that this problem is solved. Polar shapelets then pro- 
vide a natural representation for both shear and flexion operations, 
with simple mathematical forms that yield transparent physical in- 
terpretation. The complex number approach used throughout polar 
shapelets matches very conveniently with the complex ellipticity 
notation of Blandford (1991) now ubiquitous in shear literature, 
and with the complex formalism of flexion developed by Bacon et 
al. (2006). A complete software package to decompose images into 
polar shapelets is available from the shapelets web site*. 

We then exploit the inherent symmetries of polar shapelets 
to explore a comprehensive range of passive measurement meth- 
ods for both shear and flexion. To create a shear or flexion esti- 
mator, we simply need to find a combination of shapelet coeffi- 
cients that has the desired properties under each transformation. 
We generally keep the estimators as close as possible to linear in 
the image, to minimise both noise and bias in the final result. The 
shapelet methodology resembles a continuation of the KSB method 
to higher order. However, the inclusion of higher order shape infor- 
mation, and a complete parametrization of galaxy morphology, pro- 
vides several new opportunities to improve on KSB, and to remove 
its instabilities. Some of the shear and flexion estimators that we 
describe are also independent, and can be obtained simultaneously 
for each galaxy. These will provide invaluable new cross-checks for 
systematics in the data analysis, which are unique to this method, 
and can also be combined to increase the overall ratio of signal to 
noise. As we shall discuss, one of the shear estimators has already 
been proved highly successful in a blind test on simulated images 
containing an applied shear, as part of the STEP programme (Hey- 
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mans et al. 2006; Massey et al. 2007a). We defer detailed testing of 
the remainder until the next STEP cycle. 

This paper is organised as follows. In §2 we describe the 
shapelet decomposition and the action of weak gravitational lens- 
ing in shapelet space. In §3 we derive several possible weak shear 
estimators, and discuss the performance of a key estimator on the 
simulated STEP images. In §4 we derive several possible weak flex- 
ion estimators. We conclude in 85. 



2 WEAK GRAVITATIONAL LENSING IN POLAR 
SHAPELET SPACE 

We shall first describe the action of weak shears and weak flex- 
ions in polar shapelet space. This is seen as a mixing of power be- 
tween an object's various shapelet coefficients, or equivalently how 
much those coefficients change under each operation. To first order, 
a vector of shapelet coefficients is acted upon by simple matrices 
that contain small mixing components in their off-diagonal terms. 
For example, a shear takes some power from the circular (m = 0) 
shapelet coefficients and redistributes it into the elliptical (m = 2) 
shapelet coefficients, turning a circle into an ellipse. 

The effect of shear as an abstract coordinate transformation 
has already been derived in Cartesian shapelet space by Refregier 
(2003), and in polar shapelet space by Massey & Refregier (2005). 
Here, we review this shear in the physical context of weak grav- 
itational lensing. Operators to perform flexion have been derived 
in Cartesian shapelet space by Goldberg & Bacon (2005). Here, we 
translate those results into polar shapelet space, where they become 
much simpler. The flexion operators fit naturally into the complex 
notation of polar shapelets. Furthermore, the two distinct types of 
flexion identified by Bacon et al. (2006) mix distinct sets of polar 
shapelet coefficients, which can be separated elegantly. 



2.1 Polar shapelet space in the absence of lensing 

The observed image of every galaxy f(r, 6) can be decomposed 
into a sum of (complex) orthogonal 2D basis functions 
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The basis functions, which are illustrated in figure 1, are fully 
described in Massey & Refregier (2005) and Bernstein & Jarvis 
(2002). They are Laguerre polynomials in r multiplied by sines and 
cosines in 9, and a circular Gaussian of width /3. This scale size is 
chosen to match the observed size of each galaxy, and the functions 
are placed at the galaxy's centre of light. The shape of each galaxy 
can then be completely described by the array of its shapelet coeffi- 
cients fn,m- These are complex numbers, with f n ,-m = fn, m - The 
indices n and m correspond to the numbers of radial and tangential 
oscillations respectively: n can take any nonnegative integer, and 
m can take any integer between — n and n, in steps of two. The 
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Figure 1. The polar shapelet basis functions, with indices n and m that 
describe the number of radial and tangential oscillations. The functions are 
complex, but several symmetries exist to ensure that a reconstructed im- 
age is wholly real, and these have been used to condense the plot. Basis 
functions (and shapelet coefficients) with opposite signs of m are complex 
conjugate pairs. Only the real part is shown here for basis functions with 
m > and only the imaginary part for those with m < 0. The basis func- 
tions with m = are wholly real. Units of the colour scale assume that 
13 = 1. The boxes have also been enlarged into the spaces between allowed 
coefficients for clarity. 



index m will be the most significant in this paper, because coeffi- 
cients with the same value of m describe features of a galaxy with 
the same degree of rotational symmetry. 

In practice, the shapelet expansion must be truncated, and we 
typically use coefficients with n less than some maximum amount 
or, conveniently in this context, n + \m\ less than some amount. 
The latter, "diamond"-shaped truncation scheme is a cut in the to- 
tal number of oscillations, so is more consistent with arguments 
concerning information content in Fourier space, like the # m i n and 
#max of equation (24) in Refregier (2003). It is also better matched 
to the empirically observed distribution of power in shapelet space 
for typical galaxies. In figure 1, the absolute values of coefficients 
with n = 7, 8 or 9 and low |m| (which are not shown) would typ- 
ically be higher than those towards the top-right and bottom-right 
of those that are shown. For galaxy shapes this truncation scheme 
therefore improves the data compression ratio, or the accuracy of 
image recovery using a fixed number of free parameters. 

In the absence of lensing, we first assume that galaxy shapes 
are randomly oriented. This must be true for a sufficiently large and 
widely-separated ensemble of galaxies, if there is no preferred di- 
rection in the universe, and if galaxies are not intrinsically aligned. 
The unlensed ensemble of galaxies can not contain any angular in- 
formation, so must therefore have mean shapelet coefficients f nm 
that obey 



Thus, only the m = coefficients of the ensemble average are 
populated. This is the only information available about an unlensed 
galaxy ensemble. It encodes the galaxies' flux 
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and radial profile (see Massey & Refregier 2005), including their 
average size 
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and higher order shape moments like 
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as defined by Okura, Umetsu & Futamase (2006). All of these will 
be used later. 

Although the following quantities will be zero on average for 
the population, for each galaxy we can also define an unweighted 
centroid 
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the numerator of which is the /3-invariant quantity Q obtained by 
setting s = 4 and m = 3 in equations (56) and (58) of Massey & 
Refregier (2005). 

2.2 Effect of shear in shapelet space 

As a bundle of light rays from a distant galaxy passes through a 
foreground gravitational field characterized by the lensing poten- 
tial \l/(a;, y), the rays are differentially deflected, and the apparent 
shape of the galaxy is distorted (c.f. Bartelmann & Schneider 2001). 
The shape of the galaxy f(x, y) is sheared by an amount 
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Positive values of the real part, 71 , correspond to elongations 
of the galaxy along the x-axis and compressions along the y-axis. 
Positive values of the imaginary part, 72, correspond to elongations 
of the galaxy along the line y = x and compressions along the 
line y = —x. In both cases, negative values indicate the opposite. 
This complex shear notation (and an analogous form of complex 
ellipticity) is useful in weak lensing because both components are 
expected to be zero on average in the absence of a signal. In this 
case, a modulus-argument form for shear would have a zero mod- 
ulus, but no well-defined angle. The complex form also arises very 
naturally in polar shapelet space. 
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As shown in Massey & Refregier (2005), under a weak lensing 
shear S to first order, the shapelet coefficients f nm transform as 



<5 : fn,m y fn,m — fn,m (11) 

+ j | y/(n + rn) (n + rn - 2) f„- 2 ,m-2 

- \J (n - rn + 2) (n - rn + 4) f n+2 , m -2 } 
+ | \/(™ - m) (n-m- 2) /„_ 2 ,m+2 

- \/ (n + rn + 2)(n + rn + 4) / n+ 2,m+2 j 

where the asterisk denotes complex conjugation. For an intrinsi- 
cally circular galaxy, or a galaxy ensemble whose unlensed coeffi- 
cients (fnm) obey equation (3), the lensed coefficients (f' nm ) are 
left unchanged 

(/;, m ) ~ (/„,„) if m ± ±2 , (12) 
except for the \m\ —2 modes, where 
, y/njn + 2) 

Un,2> - ^ (/n-2,0 - /n+2,o) 7 , (13) 

with n = 2, 4, 6, . . .. After lensing, the galaxy has nonzero m = 
and \m\ = 2 coefficients (but no others). Figure 2 illustrates the 
action of mixing between nearby shapelet coefficients. The most 
obvious consequence is that the galaxy's unweighted ellipticity (8) 
also becomes non-zero. However, the fractional amount by which 
it changes depends upon the galaxy's radial profile. This idea will 
be explored in §3, along with other combinations of combinations 
of m = 2 coefficients. 

Note that even a pure shear to first order can change the size of 
a galaxy, if it is not intrinsically circular. But propagating series (5) 
through operation (11), and comparing the result to series (8), it is 
easy to deduce that 

S:R 2 ->R 2 ' = i? 2 (l+7£*+7*£) = i? 2 (l+27iei+272£ 2 ) .(14) 

In fact, there are (only) two different linear combinations of 
shapelet coefficients that are invariant under a first order shear: 

ri = (47T)5/3 53(/o,o + /4,o + /8,o + ...) (15) 

r 2 = (4tt)^ (/2.0 + h,o + /io,o + ...). (16) 

Furthermore, their sum is the total flux F , whose measurement is 
also independent of the choice of scale size (3. 

2.3 Effect of flexion in shapelet space 

If the shear field varies significantly across the width of an object, 
one side is distorted more than the other, and it becomes bent into 
an arclet. This effect has been dubbed "flexion". Building upon the 
work of Goldberg & Bacon (2005), we shall now describe the dis- 
tortions that arise from such gradients in the shear field, |2 . The 
calculations will remain in the weak lensing regime, in the sense 
that no terms of order 7 2 will be considered. However, flexion is 
most apparent along lines of sight close to foreground mass con- 
centrations, where the shear is also likely to be strong. The more 
rapid falloff of a flexion signal as a function of distance from fore- 
ground mass can be used to probe smaller physical scales than a 
weak shear analysis, which produces relatively non-local mass re- 
constructions. Bacon et al. (2006) demonstrate that it can be used 
to more precisely measure substructure of dark matter halos, and 
their inner profile or concentration. 



Bacon et al. (2006) pointed out that the flexion signal can be 
split into two separate (complex) terms, the first and second flexions 

^ = (jL ~ * Jy) 7 = '- 71 ' 1 + 72,2 ' ) + ^ 72,1 ~~ 7l ' 2 ' ) ^ 

^ = (Jjx + i ~&y) 7 = ^ 71,1 ^ 72 ' 2 '' + i ^ 2 ' 1 + 7l - 2 ) • 

We assume that these have the same units as 1//3 which, in the 
public code, is always expressed in terms of image pixels. Via a 
derivation analogous to that in Cartesian space by Goldberg & Ba- 
con (2005), we can determine the action of the flexion operators T 
and Q in polar shapelet space. These are much simpler than corre- 
sponding expressions in Cartesian shapelet space, because distinct 
sets of coefficients are coupled in polar shapelet space by the two 
operations, and the flexion also fits naturally into our current com- 
plex notation. 

J~ '• fn.m > fn,m = fn,m 
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where the asterisk denotes complex conjugation. Similarly, 
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These operators are illustrated graphically in figure (2). 

One crucial difference from the shear operator is that applying 
a flexion shifts the galaxy's observed centroid (7) by an amount 

R 2 

A= — {QT + ^'e + ge") , (21) 

in units of (3, with the real part corresponding to the x direction 
and the imaginary part to the y direction. The elements of expres- 
sion (21) are easily understood in terms of shapelet coefficients. A 
galaxy's centroid is constructed from its m — 1 coefficients. These 
coefficients are altered during a first flexion T if the galaxy has 
power in any m = or |mj = 2 coefficients. The m — coeffi- 
cients are never all zero, so the centroid will always shift. The cen- 
troid is altered during a second flexion Q if the galaxy has power in 
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Figure 2. The mixing of polar shapelet coefficients under weak lensing 
transformations. If a galaxy initially contains power in its f§fi coefficient, 
it will contain additional power in /4,±2 and fs,±2 after shear. After both 
types of flexion, it will contain additional power in eight shapelet coeffi- 
cients, as illustrated. The directions in which power moves between ad- 
jacent coefficients are the same for a given operator wherever there are 
non-zero coefficients across shapelet space, although the amount of mix- 
ing varies. Wherever the pattern would seem to couple coefficients that do 
not exist, the amount of mixing is zero. 



any \m\ = 2 coefficients, but the effects of its \m\ =4 coefficients 
happen to cancel out in summation (7). Therefore an object's ellip- 
ticity uniquely determines this centroid shift. No comparable shift 
was introduced during shearing, so dealing with this will present a 
new technical challenge for weak lensing measurement. 
One mapping that will be required later is 



where 



P = 



^J2( n + Wfa - 1)<> + l)(n + 3)/„ 



(22) 



(23) 



Operators (19) and (20) are useful for applying an artificial 
flexion to an unlensed galaxy (for example, during the manufac- 
ture of simulated images). However, for a practical, passive flex- 
ion measurement method, the natural location for the centre of a 
shapelet decompostion is the post-lensing (observed) centre of light 
x c , it being impossible to predict the pre-lensing sky position of the 
source. This point will be crucial in our later analysis because, for 
example, determinations of ellipticity and particularly flexion de- 
pend upon the origin of the coordinate system. To ensure that we 
account for this centroid shift, we are greatly aided by the linear de- 
pendence of operator (21) upon the coefficients that will make up 
our flexion estimators. The change in coordinate frame can be si- 
multaneously corrected for by simply incorporating an appropriate 



translation in the operator used for flexion estimation 
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where, from Massey & Refregier (2005), the translation operator is 

-^(^) ' fn,m y fn.m = fn,m (25) 
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These practical flexion operations for analysis of observed images 
effectively isolate the observable, shape-changing part of the flex- 
ion transformation by subtracting off the centroid shift. 

As described in Goldberg & Bacon (2005), for the purposes 
of constructing workable flexion estimators the ellipticity e can be 
estimated from the lensed galaxy image even though it will itself 
have changed during the lensing. The change in the centroid shift 
this represents is small, which can be seen from equation (21), and 
such changes will cancel on average due to the differing rotational 
symmetries of 7, T and Q. If deemed necessary, an estimate of 
the ellipticity corrected for locally measured shear could even be 
used, as there is nothing to prevent the galaxy shear analysis from 
being independently performed prior to any flexion analysis. These 
operators will be used to form flexion estimators from observed 
galaxy shapes in §4. 

2.4 Effect of convergence in polar shapelet space 

Convergence changes a galaxy's size and brightness. Actually mea- 
suring convergence is difficult because galaxies are intrinsically of 
very different sizes and magnitudes, and it is very hard to know 
what these quantities would have been before lensing, even statis- 
tically. (Measurements of shear and flexion are made possible by 
the statistical assumption that an unlensed population of galaxies 
would be round.) However, it is important to take account of the 
effect of convergence on these measurements, which is given by 

'^_ + dH\ 

dx 2 dy 2 ) 

Increases in apparent galaxy size potentially cause ellipticities 
to be measured in different parts of a galaxy's profile - further to- 
wards the core or out in the wings. This is compensated for by the 
adaptative choice of the shapelet scale size j3 during the shapelet 
decomposition described in Massey & Refregier (2005). Indeed, 
the operators K and S are commutative. Changes in galaxy flux, 
or the averaging of shear estimators from bright and faint galax- 
ies, can be controlled by constructing estimators that are invariant 
to object flux. This is trivially implemented for all of the estima- 
tors discussed in this paper by dividing by the flux. To first order 
in 7, this quantity is invariant under a shear. It is also the most eas- 
ily measured, zeroth-order aspect of morphology: very important 
since this appears on the denominator of shear estimators, where 
noise can translate into biases overall. 

Note that this does not mean that the issues of "reduced shear" 
(Bartelmann & Schneider 2001) or indeed "reduced flexion" (c.f. 
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Okura, Umetsu & Futamase 2006) have been solved. Pure gravi- 
tational shear or flexion are not observable in isolation. It is only 
possible to measure a degenerate combination of the shear or flex- 
ion with additional terms including the convergence. For the un- 
weighted shear estimator 7unwci g htcd, which is described in §3.4, 
the observable quantity is 7/(1 — n). However, as shown in ap- 
pendix A, this represents a limiting case that no longer holds for ar- 
bitrary weighting schemes. For convenience, the observable shear 
distortion will be labelled 7 hereafter in this paper; it should be un- 
derstood that this really refers not to the gravitational shear but to 
the reduced shear g corresponding to the estimator in question. In 
practice these reduced shears will be close to the g = 7/(1 — k) 
for the limiting unweighted case, but in appendix A we discuss how 
shapelets might be used to calculate the generalized reduced shear 
for each shear estimation method. 



2.5 Effect of convolution in polar shapelet space 

Galaxy shapes also change during convolution with a telescope's 
point spread function. In shapelet space, convolution is another 
simple matrix operation (Refregier & Bacon 2003). Deconvolu- 
tion can be performed via a matrix inversion or simultaneously 
with shapelet decomposition via a method presented in (Massey & 
Refregier 2005). We shall not further discuss the challenge of de- 
convolution in this paper, leaving it as a separable, and essentially 
solved, problem. The main effect of deconvolution is to correlate 
shapelet coefficients (since the basis functions no longer remain 
completely orthogonal after convolution). The full covariance ma- 
trix can easily be obtained during decomposition. It could, in prin- 
ciple, be used to perfect the weights on coefficients in the shear es- 
timators, although we have derived results only in the limit where 
the covariance is nearly diagonal - which is approached by basis 
functions with oscillations larger than the PSF size. 



3 SHEAR ESTIMATORS 

To measure weak shear, we would like to construct some combi- 
nation of each galaxy's observed shape components that is related 
to the shear field it has experienced. The combination can be of ar- 
bitrary complexity. For individual galaxies, the measured quantity 
will inevitably be noisy, because galaxies have their own intrin- 
sic shapes, which are changed only very slightly by weak lensing. 
However, we shall aim to construct a shear estimator 7 for which 



<7>=C> 



(27) 



when averaged over a large galaxy ensemble in the absence of 
shear; and, more importantly, 



S : 7 — > 7 + 7 



(28) 



individually. As discussed in §2.1, the first condition is easy to 
achieve by making sure that (the numerator of) 7 contains only 
shapelet coefficients with m ^ 0. The second, calibration of the 
shear estimator, ensures that the estimator is always unbiased 



Of) = 7 



(29) 



but this is notoriously difficult to satisfy (c.f. Bacon et al, 2001; 
Erben et al, 2001; Heymans et al, 2005; Massey et al, 2007a). 
Our effort will primarily be directed here. 

The easiest methodical approach towards a passive shear esti- 
mator is to first construct a "polarisation" estimator p with the same 



rotational symmetries as shear. We then need to calculate its "shear 
susceptibility" 

dpi 



/'.. _ 

so that 



S : Pi - Pi + P^lj [ + 0(7 2 )] 



(30) 



(31) 



The shear susceptibility can usefully be thought of as two complex 
numbers; one for each component of shear. However, it is more 
commonly expressed as a real, 2x2 tensor and, for the sake of fa- 
miliarity, we shall adopt that notation here. Its diagonal (real) terms 
describe the amount by which the polarisation will change under a 
shear. The off-diagonal (imaginary) terms describe a peculiar mix- 
ing by which a shear in one direction can affect the polarisation in a 
direction at 45°. This is introduced by complex galaxy morpholo- 
gies when a galaxy's isophotes are not concentric. 
We can then construct a shear estimator 



7i = (P^i/Pi ■ 

to make sure that indeed 

Hi) = mr'pj + (p^r'p^i) 

= <(^)- 1 ^> + <7i> 
= 7; i 



(32) 

(33) 
(34) 
(35) 



where the random intrinsic ellipticities of galaxies ensure that the 
first term vanishes, and thus condition (29) is satisfied. 

However, we immediately encounter four difficulties with 
shear susceptibilities that account for most of the problems in the 
current generation of shear measurement methods: 

• P 7 is noisy. It is usually constructed from an object's higher 
order shape moments, which are even harder to measure than the 
polarisation. Since this appears on the denominator, it dramatically 
increases the scatter of the shear estimator: any ratio of quantities 
with Gaussian errors produces the extended wings of a Cauchy dis- 
tribution (as seen for a KSB analysis in figure 2 of Massey et al. 
2004), whose moments like 07 do not even converge. 

• P 7 is a tensor. The matrix inversion in equation (32) is unsta- 
ble, except for circularly symmetric galaxies, or an unlensed popu- 
lation ensemble, in which case the off-diagonal elements are always 
zero. In all other cases, shearing in one direction mixes ellipticity 
from all other directions, and this must be unmixed. 

• P 7 is required pre-shear. Each galaxy is observable only af- 
ter it has been lensed. Unfortunately, the shear susceptibility factor 
may change during shear, to first order in 7 for most galaxies, and 
to second order for even circularly-symmetric ones. 

• The P 7 formalism ignores terms of second order in shear. 
This omission may bias shear measurements at the sub-percent 
level of precision, and introduce non-linearities that depend upon 
an object's intrinsic ellipticity and \m\ = 4 shapelet coefficients. 

A frequently adopted solution to the first three difficulties is 
to average P 1 from a set of intrinsically similar galaxies, or to fit a 
value from a large galaxy ensemble as a function of other observ- 
ables. This approach ought to find a suitable, statistical value for 
all galaxies. It diagonalises the shear susceptibility; reduces noise; 
and, if the population is so large that it contains effectively no co- 
herent shear signal, satisfies the requirement for the pre-shear mea- 
surement. The fourth difficulty is particularly troublesome because 
an object's measured ellipticity is degenerate with the shear - but 
may also be resolvable in averages over a large population of galax- 
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ies chosen without shear-dependent biases. Unfortunately, averag- 
ing over any large population of galaxies is inelegant, in the sense 
that shear estimators for individual galaxies are no longer self- 
contained. It also introduces new problems: the main issue being 
the practical identification of a set of intrinsically similar galaxies. 
Most observable properties of a galaxy do change during a shear, 
and grouping galaxies by these leads to "Kaiser flow" (Kaiser, 
2000). The common challenge facing all modern shear measure- 
ment methods is to either understand Kaiser flow statistically, or 
to control shear susceptibility and thus avoid it. In appendix B, we 
show how measurements with one polarisation estimator can be av- 
eraged to avoid Kaiser flow, and maximise the weak lensing signal. 

For the rest of this section, we shall construct progressively 
more elaborate polarisation estimators that ameliorate the four dif- 
ficulties. We begin with simple polarisations that are compactly 
represented in polar shapelet space. These still suffer from all four 
difficulties. We then gradually exploit the symmetries of shapelets 
to add more complex features. The process is helped by the conve- 
nient shapelet notation, although the expressions do become more 
complicated. Which of these advanced shear estimators is most ap- 
propriate to a given data set will depend on the desired applica- 
tion, the image quality (for example, whether it was taken from the 
ground or in space), and the number of shapelet coefficients avail- 
able for each galaxy. 



sion, i.e. whether to normalise by flux or not, will also have to be 
made for all of the following shear estimators. 

3.2 Order-by-order shapelet shear estimator 

A successful shapelet decomposition contains all of the available 
information about a galaxy's shape, and more information can be 
extracted than that available with previous shear estimators. Since 
all of the \m\ — 2 shapelet basis functions have the same rota- 
tional symmetries, each of the corresponding shapelet coefficients 
can be used to form independent (except for the covariance be- 
tween shapelet coefficients after deconvolution) polarisation esti- 
mators p — f n ,2- These have shear susceptibilities 

(PZ)n = \{Vn(n + 2) (/„- 2 ,o - /n+a,o) (40) 
+ v/(n-4)(n-2) Re {/„- 2 , 4 } 

- ^(n + 4)(n + 6) Rc{/„ +2 , 4 }} 

(PZ)22 = i{V»»(n + 2) (/n-2,0 - /n+2, ) (41) 

- ^(n-4)(n-2) Rc{/„_ 2 , 4 } 
+ ^(n + 4)(n + 6) Re{/ n+2 , 4 }} 



3.1 Gaussian-weighted quadrupole moment 

We shall start with the simplest possible combination of shapelet 
coefficients that can be used to build a polarisation estimator. Re- 
call that the first shapelet coefficients to be affected by a shear are 
those with \m\ = 2. Like shear, these rotate as e -2 *, and they are 
therefore suitable for our purposes. The simplest possible polarisa- 
tion estimator is simply the first shapelet coefficient with m = 2, 
i.e. p — / 2j2 . This has shear susceptibility 



Pii = (/o,o - U,o)/V2 - V3Re {/ 4 , 4 } 
P22 = (/o,o - /4,o)A/2 + \/3Re {/ 4 , 4 } 
P12 = P21 = -\/3Im{/ 4 , 4 } . 



(36) 
(37) 
(38) 



In images from the Hubble Space Telescope COSMOS survey 
(Scoville et al, 2007), for example, (|/ 4 , 4 |//o,o) « 0.079, which 
is not entirely negligible at the desired level of precision. By av- 
eraging the components of P 1 from a sufficiently large population 
of observed galaxies, or fitting them as a function of other observ- 
ables like galaxy size and magnitude, we can explicitly force the 
mean m — 4 coefficients to be zero, and ensure that the measured 
m — coefficients are statistically corrected before shear. With 
this simplification, the shear susceptibility factor can then be triv- 
ially inverted, and we arrive at the shear estimator 



7G. 



aussian 



(/o,o - / 4 ,o) 



(39) 



This recovers the methods of RRG (excluding the smear cor- 
rection) and Refregier & Bacon (2003), casting them into the more 
succinct framework of polar shapelets. It recovers the P sh compo- 
nent of KSB up to the normalisation of the polarisation estimator. 
To avoid biases and instability at low signal to noise, we have cho- 
sen to keep the polarisation and shear susceptibility linear in the 
image brightness. As a result however, both quantities vary widely 
in the full galaxy ensemble which typically encompass large ranges 
of flux and sizes, increasing the rate of Kaiser flow. A similar deci- 



{P2h 



(PS) 



= i| N /( n -4)(n-2) Im{/„- 2 , 4 } (42) 
- V(n + 4)(n + 6) Im{/„ +2 , 4 }} , 

which reduce to 



u-f _ Vn(n + 2) 

— -. U n-2,0 - J n+2,0) 



(43) 



when averaged over an ensemble of galaxies as before. Thus, for 
each even order n available in a shapelet decomposition, we can 
construct one independent, unbiased shear estimator 



In = 



f'n,2 



, for n = 2, 4, 6, . . . (44) 



\Jn(n + 2) (/n-2,0 — j "n+2,0) 

As before, these estimators are by construction unbiased, when av- 
eraged over the galaxy population. 

One way to use these additional estimators is to diagnose prob- 
lems in the measurement. Because we obtain multiple shear mea- 
surements for each galaxy during a single PSF deconvolution, their 
agreement provides a strong new test of systematics. If a pure shear 
signal is being successfully measured, all of the estimators from a 
given galaxy should average to the same value. However, if resid- 
ual PSF effects are polluting the signal, the separate estimators will 
disagree. A weak lensing pipeline must be highly robust to pass 
such stringent tests, and they will provide a unique discriminatory 
power in future analyses. 

Alternatively, the separate estimators can be linearly com- 
bined, with arbitrary weightings 



(45) 



where the summation runs only over even indices n, for only those 
coefficients exist. In this case 



(46) 
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The weights w n can be carefully constructed to optimise the 
signal-to-noise of the shear measurement (such as inverse variance 
weighting, as suggested by Refregier & Bacon 2003) or to remove 
systematic biases plaguing the particular data set. For the rest of 
this section, we shall explore various options for this weight func- 
tion. By staying linear in shapelet coefficients during this process, 
the polarisation and susceptibility also stay linear in the image, thus 
preserving a Gaussian-like distribution of estimators. In real space, 
changing the weights w n is equivalent to changing the weight func- 
tion used for the polarisation estimator. 



3.3 Using galaxies' radial profiles to reduce a 7 

A galaxy's observed m = 2 coefficients consist of intrinsic ellip- 
ticity, shear-induced ellipticity, and noise. For an individual galaxy, 
there is no way to tell what fraction of each is intrinsic, and what 
fraction is the signal. However, a shapelet decomposition contains 
a great deal more information about a galaxy's morphology that has 
not yet been tapped. In particular, it is the galaxy's intrinsic radial 
profile (m = coefficients) that contribute most to any change in 
observed ellipticity during a shear. Since the m = coefficients 
are typically much larger than any others, they are also, fraction- 
ally, the least changed themselves under a small shear. We shall 
therefore approximate the unlensed radial profile as the observed, 
measured radial profile. We can then work out the "radial profile" of 
m — 2 coefficients that could possibly have been induced by lens- 
ing. Any component of the intrinsic ellipticity that does not have 
the appropriate "radial profile" cannot possibly have been induced 
by lensing and, for our purposes, can be ignored. Thus we reduce 
the contamination of intrinsic galaxy ellipticity in our shear estima- 
tors, to only include components of intrinsic ellipticity that happen 
to have the right profile. 

We determine the required weights w„ by applying a unit 
shear to the rotationally-invariant part of a model, and find 



7profilc 



L S V n ( n + 2) (/n-2,0 ~ /n+2,p) fn,2 

(£n(n + 2) (/„_ 2 ,o - / n+2 ,o) 2 ) 



(47) 



where one factor in the denominator comes from the shear suscep- 
tibility factor and one from the weighted average. Of course, we 
have not taken measures to eliminate the \m\ = 4 and off-diagonal 
terms in the shear susceptibility factor. The shear susceptibility will 
therefore need to be fitted from a galaxy population as a function of 
size, magnitude and possibly radial profile. Several shapelet-based 
parameters to span morphology variation are suggested in §7 of 
Masssey & Refregier (2005). 



3.4 Diagonal shear susceptibility 

One of the difficulties with general shear estimators, as described 
at the start of §3, is that they require the inversion of a (noisy) 
shear susceptibility tensor (46). This inversion is often unstable, 
and various implementations have chosen to either ignore the off- 
diagonal elements, or average over a large population of galaxies 
so that they disappear. The problem could be solved more easily if 
the shear susceptibility were explicitly a simple scalar (times the 
identity matrix) for each galaxy. Indeed, it is possible to weight the 
various orders of 7„ in such a way that the off-diagonal terms in 
their combined susceptibility tensor from successive orders cancel 
each other. The off-diagonal terms, and the differences between the 
on-diagonal terms, involve \m\ = 4 coefficients that are introduced 
by the 7* terms in equation (11). With these removed, the shear 



susceptibility (46) will be diagonal and only involve terms with 
m = 0. This can be trivially inverted. 

A simple calculation to obtain the desired w n yields 



p = \M™ + 2) fn2 



(48) 



where the prefactor has been added to reproduce familiar quanti- 
ties. In fact, p — FR 2 e, a version of the (radially) unweighted 
ellipticity without size (or flux) normalisation. This can not nor- 
mally be calculated from images because background noise makes 
the real-space integrals diverge. A shapelet decomposition removes 
noise by acting as a prior on the permitted physical properties of a 
galaxy shape. This polarisation has shear susceptibility 



P 7 = 16ttV jr VnTT Uo = 2FR 2 



(49) 



where the right hand side refers to quantities measured before 
shearing. The susceptibility is the size and magnitude of galax- 
ies, in a curious contrast to the previous shear susceptibilities that 
needed to be ensemble averaged and fitted as a function of those 
observables. Furthermore, as shown in equation (14), to first order 
in 7, the size R 2 changes under a shear in a way that affects the 
overall shear estimator 



p 



p + 2FR 2 ~f 



(50) 



' 2FR 2 2F'R 2 ' 2Fi? 2 (l + 7£* + 7*6 
Ensemble averaging, and expanding to first order in 7, we recover 

(e 2 ) 



P 

2F'R 2 ' 



\2FR 2 



(51) 



(2) 

with the same "shear responsivity" factor of 1 — that appears 
in equation (A14). Thus we obtain an unbiased shear estimator 



7unwcightcd 



£Vn(n + 2)/„ 



(2 - <£ 2 >) E \ftn + Tjfn0 



(52) 



that is written in terms of observable quantities alone, and requires 
minimal averaging of shapelet coefficients from a population of 
galaxies. 

This particular shear estimator emerged as one of the most 
successful shear measurement methods during blind tests as part 
of the STEP programme (Massey et al. 2007a). This programme 
constructed simulated images that exhibit all the statistical prop- 
erties of real astronomical images, but contain a known shear sig- 
nal. While the measurement was performed (by JB), these input 
shears were kept hidden. They were then revealed publicly after all 
the pipelines had been run. Figure 3 shows the impressive perfor- 
mance of our Tunwcightod shear estimator for STEP2 image set A, 
which was specifically designed to mimic deep Suprime-Cam im- 
ages from the Subaru telescope (Miyazaki et al., 2002). A linear 
fit to these results shows a shear calibration factor (multiplicative 
measurement bias) of m = 0.023 ± 0.029 for the real component 
of shear, and m = 0.053 ± 0.029 for the imaginary component. 
There is no significant residual shear offset (additive measurement 
bias), with the fitted values being c = (—6.8 ± 6.5) x 10~ 4 for the 
real component of shear and c = (1.3 ±6.6) x 10~ 4 for the imagi- 
nary component. This estimator demonstrated the best performance 
throughout the STEP2 project and, as suggested by that analysis, 
we shall reduce our error bars before the next round by incorporat- 
ing a galaxy weighting scheme into our weak lensing pipeline. 

One additional benefit of this estimator is that both the polar- 
isation and the shear susceptibility are independent of the shapelet 
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Figure 3. Performance of the shear estimator with an explicitly diagonal 
shear susceptibility tensor, on simulated images containing a known true 
shear. Points show both components of the mean shear, measured in 7' X 7' 
patches of sky where the input shear was held constant. For a perfect shear 
estimator, these would lie on the solid line. The dashed line is a linear fit to 
deviations from it. 

scale scale (3. Although the ensemble average of any shear esti- 
mators (7„) should always be independent of f3, in general, indi- 
vidual estimators 7„ may not be. But in the current case, once a 
shapelet series has converged, F, R 2 and e combine coefficients 
in such as way as to not depend upon the choice of j3 (Massey & 
Refregier, 2005). This result is non-trivial: in our image decomposi- 
tion pipeline, we choose f3 to optimise the image reconstruction and 
stabilise the PSF deconvolution. However, this is only one possible 
goal. In IMCAT implementations of the KSB method, the equivalent 
of /3 is instead chosen to maximise the signal-to-noise for detection 
of the object. In SEXTRACTOR implementations of KSB, a scaling 
of SEXTRACTOR parameters is used. In all cases, the choice of f3 
will also be affected by any applied shear that changes a galaxy's 
apparent size. Whether this change is negligible depends on the 
specific implementation of the algorithm to determine f3 and must 
be tested experimentally. In this case, it is reassuring to note that 
this effect is formally absent, modulo the convergence of the series. 
Since the STEP tests, we have also derived the full "Kaiser 



flow" behaviour of this estimator. This has a particularly interesting 
form, and is discussed fully in appendix B. 



3.5 Shear-invariant shear susceptibility 

The shear susceptibility can be further simplified, by continuing 
to add terms to the polarisation estimator. So far, we have used 
the |m| = 2 coefficients, but no others. In fact, all shapelet ba- 
sis functions with \m\ = {2, 6, 10, 14, . . .} contain the rotational 
symmetries of shear; the higher order functions just contain addi- 
tional symmetries as well. They can not be used as shear estimators 
by themselves but, because of this, it is possible to add them to the 
m | —2 coefficients. The resulting polarisation estimator will stay 
linear in shapelet coefficients, linear in image flux, and keep all of 
the desired symmetries. Here we shall demonstrate how higher or- 
der shapelet coefficients can be used to "sweep" terms in the shear 
susceptibility down to to = 0, and construct a polarisation 

oo oo 

P = ^ W n,mf n ,m (53) 

m=2,6,10,... n=2,4,6,... 

with any desired susceptibility factor. 

We begin with the /2,2 polar shapelet coefficient. As shown 
in §3.1, this has a shear susceptibility involving /o,o, /4,o and /4,4. 
The weight on the /2,2 coefficient in p can be used to create any 
desired factor in front of the / ,o term in P 1 . In §3.4, we added 
coefficient /e,2 in such a way to cancel out the /4,4 coefficient in 
the shear susceptibility. Instead, we shall now add an amount of 
/6,2 (then /io,2, /i4,2, etc) to shape the susceptibility's /4,o (and 
/s,o, /12.0) terms in any way. Assuming extrapolation to infinite n, 
we can thus construct a polarisation estimator with arbitrary m — 
terms in its shear susceptibility tensor. However, it will also contain 
non-zero \m\ = 4 terms and off-diagonal elements. 

Now consider the fe t e polar shapelet coefficient. This has a 
shear susceptibility involving /4,4, fs,4 and /s,8 coefficients. This 
can be added to the polarisation with a weight arranged so that 
the three /4,4 terms in the shear susceptibility now cancel out. The 
/io,6 coefficient can then be added so that the four /s,4 terms can- 
cel, and so on. This leaves 'dangling' terms in the shear suscepti- 
bility with \m\ = 8 (and due to series truncation, in practice, high 
n). Successive additions of \m\ = {10, 14, 18 . . .} terms to the po- 
larisation can push these terms to higher and higher \m\. Since the 
magnitude of shapelet coefficients for real galaxies typically fall off 
rapidly with n and \m\, the contribution of any remaining dangling 
terms due to series truncation decreases during this process. 

A second, interwoven combination of shapelet coefficients 
starting with /4,2, /s,2 and /s,e can also be constructed, to make 
a completely separate polarisation whose shear susceptibility in- 
volves arbitrary contributions of only /2,o, fe,o, ■ ■ ■ coefficients. 

We are now free to decide the most suitable form for the shear 
susceptibility, and can construct any appropriate polarisation esti- 
mator. It would be easiest to satisfy requirement (28) if P 1 did not 
change under shear. We could then use the observed (post-shear) 
value for each individual galaxy. The only two quantities (15) and 
(16) like this can clearly be constructed from these two interwoven 
combinations of shapelet coefficients. We shall construct the polar- 
isation estimator with shear susceptibility T\ + T2 = F, where F 
is the flux. Clearly, any shear estimator must eventually be nor- 
malised so that the same shape is calculated for two (otherwise 
identical) galaxies of different brightnesses. The flux is the sim- 
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Table 1. Coefficient weights for the real part of the 7 s hcar-invariant shear 
estimator (equation 54). The imaginary part of this shear estimator involves 
complex conjugates of terms with every other value of m, but is most con- 
veniently found by instead rotating the galaxies by 45° then calculating the 
real part a second time. 



2 
4 
6 
8 
10 
12 
14 
16 
18 




10 
12 
14 
10 
12 
14 
14 



6 
6 
6 
6 
6 
10 
10 
10 
14 



8/(3^20) 
4/(3^35) 
8/(3%/35) 
8/(3a/105) 
8/(3-v/42) 
16/(15^) 
16/(15a/77) 
96/(15a/462) 
64/(7v / 858) 



plest, and most robustly measured normalisation, and the obvious 
choice to put on the denominator; i.e. 



7shcar— invariant 



'n,m J n,m 



(54) 



Calculating w n , m for this case is an elementary but tedious 
procedure using the iterative method described above. The first 
terms are listed in table 1 . Note that although each term in the se- 
ries is well-defined, a real-space calculation in Appendix C demon- 
strates that the summation approaches a solution that is not contin- 
uously differentiable. 



3.6 Convergence issues 

Unweighted ellipticities (and higher-order properties) can not usu- 
ally be measured from real images, due to the presence of noise 
that makes the integrals diverge. However, this becomes possi- 
ble with shapelets (or IM2SHAPE or GalFit) because a noise- 
free (and unpixellated and deconvolved) model of the galaxy is re- 
constructed. For different weight functions, the question becomes 
one of the convergence of galaxies' radial profiles to zero at large 
radii, and the suitable truncation of its measured shapelet series. A 
galaxy with an exponential profile has polar shapelet coefficients 
f n fi oc n~ -° for a well chosen 13: for a sufficiently high all 
of the shear estimators do converge. For a further discussion of the 
convergence of shapelet series, see Massey & Refregier (2005). 

Furthermore, although polar shapelet coefficients can be mea- 
sured to infinite n using the linear overlap method (see Massey & 
Refregier 2005) for idealised data, this is not the case using the now 
standard x fitting method (see also Kuijken 2006; Nakajima & 
Bernstein 2006), or in the presence of pixellisation or a PSF (Berry, 
Hobson & Withington 2004). We therefore need to ensure that suf- 
ficient coefficients can be measured, particularly for the elaborate 
polarisation estimators that converge more slowly. They may there- 
fore require galaxies of slightly higher signal to noise and n ma x. 
This can be achieved by raising the magnitude or size cut in a lens- 
ing catalogue, although the extent to which this is necessary has to 
be determined by experiment. However, the more elaborate polar- 
isation estimators have correspondingly simpler shear susceptibili- 
ties, which converge faster. As was the case for noise, the minimi- 
sation of truncation errors is particularly important in the denomi- 
nator, and this may in fact prove to be the deciding factor. 



3.7 Active shear estimators with polar shapelets 

Bernstein & Jarvis (2002), Kuijken (2006) and Nakajima & Bern- 
stein (2006) suggest a rather different philosophy for constructing 
shear estimators. Instead of measuring an observed polarisation, 
and calculating how that would have changed during shear, Bern- 
stein & Jarvis (2002) shear objects (or their coordinate system) un- 
til they appear circular. Kuijken (2006) assumes that objects were 
intrinsically circular (i.e. models with f n ,m — V m 7^ 0), then 
shears them until they most closely resemble the data. This makes 
for a simpler calculation because the shape of each object, includ- 
ing its higher order moments, is known well before the operation. 
Forward shearing of an image is also useful, because it can be per- 
formed to arbitrarily high order in 7, addressing the fourth concern 
in §3. In real space, the sheared ellipticity of the shapelet basis func- 
tions in Bernstein & Jarvis (2002) can be chosen from a continuous 
range; in shapelet space, operation (11) can be exponentiated to in- 
clude higher order terms. However, a shear susceptibility factor is 
still needed (the calibration factor TZ in equation (5.33) of Bernstein 
& Jarvis (2002) has the same origin as that in our equation (50)). 
This too must be fitted or interpolated from a galaxy population as 
a function of other observables, so offers no advantage over a shear 
susceptibility. 

At first sight, this shear measurement philosphy seems more 
efficient than the passive, moment-based shear estimators consid- 
ered in the rest of this paper. Active methods do not require the full 
model of each deconvolved galaxy shape, but extract only the re- 
quired quantity 7. However, the decadence of obtaining a full shape 
reconstruction can provide extra information that is invaluable. For 
example, checking that the model's residual image is consistent 
with noise can indicate potential problems, and which shear esti- 
mates to trust, in a way that is not possible if ony a single number 
is obtained for each galaxy. In principle, it is possible to expand the 
definition of "circularity" in Bernstein & Jarvis (2002) to involve 
different shapelet coefficients, but this does not generate extra in- 
formation that is necessarily useful for systematic cross-checks. In 
that case, each fit would require a separate non-linear iteration to 
find the best-fitting parameters x c , /3, n max and 7, and therefore 
could be subject to independent biases. Since altering the shear es- 
timator could change any systematic influences in this method, it 
would be difficult to interpret any variation between the estimators. 
Instead fitting a model that is simultaneously capable of captur- 
ing all the shape information also makes the PSF deconvolution 
more robust and intuitive than methods that use a model represent- 
ing only the best-fit elliptical profile of a complex galaxy shape. 

In our experiments with elliptical shapelet basis functions, we 
have confirmed that the choice of that ellipticity is the most un- 
stable part of the iteration, particularly for objects at low signal 
to noise. We have had one idea for a different truncation scheme 
with highly elliptical basis sets. A problem arises when fluctuations 
along the minor axis become smaller than the PSF or a pixel and 
therefore non-orthogonal. Simply decreasing n max (Nakajima & 
Bernstein 2006) also shortens the reach of the basis funtions along 
the major axis, and therefore could potentially bias a shear mea- 
surement. However, it is possible to first rotate the basis functions 
so that 8 = lies along the major axis of the ellipse, then truncate 
the basis functions at different values of m and 712, the Cartesian 
shapelet indices in the x— and y— directions. If the newly-truncated 
coefficients are kept, but initially set to zero, operation (37) from 
Massey & Refregier (2005) can be used to recover the coefficients 
that would have been obtained from an unrotated set of elliptical 
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basis functions and thus continue the Bernstein & Jarvis (2002) 
method. 

4 FLEXION ESTIMATORS 

4.1 Gaussian-weighted flexion estimators 

We shall now try to develop estimators for the weak lensing flexion, 
making use of our experience with the shear estimators, and draw- 
ing tight analogies. The simplest passive simple flexion estimator 
can be formed from a similar approach to that taken when con- 
structing our KSB-like shear estimator ^Gaussian- For that, we con- 
sidered the Gaussian weighted quadrupole moments, which were 
the first perturbed under a shear. In this case, the shapelet coeffi- 
cients primarily affected by first and second flexions, transform as 



Ft : /i,i — » /i,i 



+ 



F(3 



T*(3 



6 1- 



/o,o + 6^-/2,0(55) 



- 6/ 4 ,o-5V2e* |s7a,aj 



R 2 

Se-sa (/o,o - h,o) 



+ +6^-) /a,2, -3^6/4,2 J , 



The flexion susceptibility factors are real, 2x2 tensors, and can be 
calculated using equations (19) and (20). The need to subtract away 
an estimate of the shift in the galaxy centroid due to the flexion it- 
self, expressed by equations (24) and (25), necessarily complicates 
these expressions. However, this then describes the measurable ef- 
fect of flexion on galaxy images. The first flexion susceptibility for 
general m = 1 polar shapelet coefficients is 



(Pf )n + i(Pnhi = 

( 3-vM + 1 [(n - l)(/„_ 3 ,o - /n+1,0) + 
16^2 \ ( n + 3)(/ TO _ li0 _/ n+ 3 i0 )] 



+ 3y/(n-3){n- l)(n + l) /„_ 3 , 2 
+ (3n + ll)Vn=T /n-1,2 
- (3ra - 5)Vn + 3 /n+1,2 
- 3x/(n + l)(n + 3)(n + 5) /„ +3 ,2 
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.R 2 



+ 2-^(6 + 5e r )(V^+3f n+1 ,2 - V^T/„_i, 2 ) 



(61) 
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Therefore, the combinations 

4/3 / M 



3 <(/3 2 - P 2 )/ ,o + P 2 / 2 ,o - /3 2 / 4 ,o> 



and 



4^6 



h 



(57) 



(58) 



3/3 (/ ,o + /2,o — fi,o — fe,o) 

can be used as simple flexion estimators. 

Note that the es in equations (55) and (56) refer to the pre- 
lensing ellipticity, and really do cancel out when averaging over 
a population of galaxies, even in the presence of a shear field. 
Changes in R 2 due to flexion do not bias (R 2 ) to first order either, 
as these cancel when averaged over a population of galaxies. 

4.2 Order-by-order shapelet flexion estimators 

For the small and faint galaxies that make up the majority of a weak 
lensing survey, it will be difficult to measure polar shapelet coeffi- 
cients beyond the n — 6 terms needed to estimate Q as described 
above. However, for those galaxies whose higher order shapes can 
be measured, it is possible to generalise the flexion estimators. 

The terms in curly brackets in equations (55) and (56) effec- 
tively describe a flexion susceptibility factor, which we introduce 
by analogy with the shear susceptibility factor (31). We shall then 
be able to form flexion estimators 



F n — {(Pn /n,l 



and 



fn,n 



(59) 
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The second flexion susceptibility for m = 3 coefficients is 



(P n 6 )u+i(P n e ) 2 i = 



P ^(n- l)(n + l)(n + 3) x 
16^2 | (/n-3,0 + /n-1,0 - /n+1,0 - /n+3,0) 

+ y/(n - 7)(n - 5)(n - 3) /n-3,6 
+ v/(^-5)(n-3)(n + 5) /n-1,6 

- ^(n-3)(n + 5)(n + 7) /„+i, 6 

- v/(n + 5)(n + 7)(n + 9) /n+3,6 
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(64) 



P y/{n-l)(n + l)(n + 3) X 
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In all four cases, the last six terms are complex, and the final two 
emerge from the shift in an object's apparent centroid during flex- 
ion. 

We shall now consider options by which m = 1 and m — 3 
coefficients of different orders n can be combined. We search for 
sophisticated combinations that produce flexion estimators with 
useful properties, analogous to those already created for shear esti- 
mators. 



4.3 Using galaxies' radial profiles to improve flexion 
estimators 

Exactly as was done for shear estimators in §3.3, it is possible to use 
knowledge of a galaxy's radial profile to restrict which component 
of its | m\ = 1 or \m\ = 3 polar shapelet coefficients could have 
been induced by flexion. Via a parallel derivation, we obtain flexion 
estimators 



profile — 



16\/g J2 W nfn,l 

3/3 (£K+i)) 



where 

w n = Vn + f(n - l)(/„_3,o + /n+i,o ) 
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where 
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(67) 



Wn = \/(n - 3)(n - l)(n + 1) x (68) 

(/n-3,0 + fn-1,0 — /n+1,0 — /n+3,0) • 

These are guaranteed to converge for a typical galaxy if suffi- 
cient terms are available in its shapelet series. The estimator for the 
second flexion in particular should provide a very clean measure- 
ment with minimal noise. 



4.4 Diagonal flexion susceptibility 

It might also be hoped that successive off-diagonal terms in the 
flexion susceptibility matrices could be made to cancel via a suit- 
able weighting scheme w„, as was possible for shear in §3.4. Un- 
fortunately, due to the presence of the centroid-shifting correction 



so necessary for reliable flexion estimators, this is difficult; espe- 
cially for the first flexion. 

For the second flexion we can come tantalisingly close, and 
indeed if we only consider the pure Q transformation of equation 
(20), the weighting scheme w n — y/(n — l)(n + l)(n + 3) can 
be used to form a second flexion estimator 



diagonal 



2V2 E v/(n-l)(n + l)(n + 3) f„, 3 



3/311 



E(n 2 + 2n + 2) / n , 



(69) 



This is none other than the quantity 45/3, as developed for 
HOLICS by Okura, Umetsu & Futamase (2006), except for the 
additional "flexion responsivity" factor TZ. This arises because the 
denominator changes during flexion (see equation 22), biasing the 
overall estimator by an amount 1 — -^p in a completely analogous 
fashion to the shear responsivity factor calculated in §3.4. Also, 
the inclusion of terms from the flexion-induced centroid shift (24) 
results in off-diagonal elements in P e that cannot all be removed 
through any combination of the m = 3 coefficients. 

In the case of the first flexion, the prospects are worse: even if 
we could omit the T part of a practical flexion operator (which, 
for T we most certainly can not), a w n capable of cancelling 
the off-diagonal terms in the susceptibility matrix has yet to be 
found by the authors. The complication arises from the mixing 
of power between Am, An = ±1 coefficients in the first flex- 
ion operation (19). Like a centroid shift, flexion causes power to 
leak between adjacent shapelet coefficients (c.f. figure 2). However, 
whereas the centroid shift involves only the single ladder-operator 
transformations al, a\, a r and a; (see Refregier 2003), flexion al- 
ways acts via combinations of three of these ladder operations, tak- 
ing three steps but doubling back to move only one overall. Since 
at does not commute with a r , nor a\ with ai, each Am, An = ±1 
term in equation (19) is in fact a combination of five separate con- 
tributions, each of which representing a different, independent path 
between the coefficients. Worst of all, each path contributes a dif- 
fering, n-dependant proportion of the overall transformation. This 
added level of complexity for the first flexion transformation there- 
fore precludes any estimator of first flexion with vanishing off- 
diagonal terms in the susceptibility matrix. 

The /3-invariant quantity obtained by setting s = — 1 and m = 
1 in equations (56) and (58) of Massey & Refregier (2005) could be 
used to measure first flexion. Unfortunately, this quantity does not 
appear to have any other properties that are particularly interesting 
in the context of weak lensing. 



4.5 Active flexion estimators with polar shapelets 

In a similar way to the active shear estimators, it is also conve- 
nient to use a shapelet representation when distorting a circular ob- 
ject (or possibly an object with both circular m = and elliptical 
\m\ = 2 components) by applying flexion until it matches the ob- 
served shape. Goldberg & Bacon (2005) suggested a prescription 
in Cartesian shapelets, which has been implemented by Goldberg 
& Leonard (2007), to fit the odd shapelet coefficients by perturbing 
the even ones. This results in a simple, x 2 minimisation via ma- 
trix inversion. However, their approach is not perfectly clean. The 
even Cartesian shapelet coefficients mix a great deal of structure 
beyond the circularly symmetric and elliptical components. These 
are isolated using polar shapelets and, furthermore, so are the first 
and second flexion signals. By using polar shapelets, it is possible 
to fit T and Q independently, from the \m\ — 1 and mj = 3 po- 
lar shapelet coefficients. Since the flexion signal is spread evenly 
across fewer polar shapelet coefficients than Cartesian ones, but 
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noise from an uncorrelated sky background is equal in all shapelet 
coefficients, using fewer coefficients will result in a cleaner fit, with 
improved signal to noise. 



5 CONCLUSIONS 

We have described the mechanics of weak gravitational lensing 
in terms of "polar shapelets" (Refregier 2003; Massey & Re- 
fregier 2005). This is a set of basis functions that can be used to 
parametrize images, in a way that appears convenient for both weak 
shear and flexion measeurement. The symmetries of polar shapelets 
are well-matched to those of lensing. For example, the complex no- 
tation of polar shapelet coefficients, where their modulus represents 
the amount of power, and their phase represents their orientation, 
mirrors that commonly used in the literature to define a complex el- 
lipticity. In addition, polar shapelets concisely encapsulate the ideas 
of weak flexion that have been recently developed. 

The symmetries inherent to the polar shapelet formalism pro- 
vide useful insight into the parallel symmetries of weak lensing 
distortions. We have exploited this relation to construct estimators 
that are able to simultaneously extract both the weak shear and flex- 
ion signal from the observed shapes of distant galaxies. We attempt 
to bypass some of the limitations of KSB and other shear measure- 
ment methods that were reviewed in STEP2. Adopting the clas- 
sification scheme from that programme, we briefly discussed the 
recasting of alternative, "active" shear and flexion estimators into 
the polar shapelet formalism, and more comprehensively explored 
the options available for "passive" shear and flexion estimators. 

The Gaussian-weighted shear estimator TGaussian recovers old 
methods like KSB and RRG, but in a compact complex notation. 
The unweighted shear estimator ^unweighted takes advantage of the 
noise-free shapelet reconstructions to diagonalise the shear suscep- 
tibility tensor into a scalar quantity. This particular quantity hap- 
pens to be easily derivable in real space as well. The simplifica- 
tion of the shear susceptibility is completed with 7 s hcar- invariant- 
With this, the shear susceptibility is simply the object's flux: a ro- 
bustly measured quantity, and one that does not change to first or- 
der during shear. The growing complexity of these shear estima- 
tors solves progressively more of the four issues highlighted with 
previous-generation shear measurement methods. However, they 
also converge more slowly, and require high-n coefficients to be 
available. The later estimators may consequently be accessible only 
on galaxy images with higher signal-to-noise. The best estimator to 
use (which may not even be the same for an entire population) may 
therefore depend on the flux of an object, or the nature of residual 
problems found in any particular dataset. One final, particularly in- 
teresting alternative option is the estimator 7p r0 nie that can reduce 
the rms shear noise due to intrinsic galaxy ellipticities, by explot- 
ing additional information about each galaxy's radial profile. Anal- 
ogous estimators for flexion also exist for most of these options. 

Interestingly, our method permits several independent shear 
estimators and several flexion estimators to be obtained simultane- 
ously for each galaxy. Because we calculate them following a single 
PSF deconvolution or non-linear iteration, their agreement (or oth- 
erwise) will provide a stringent new test on the PSF modelling and 
for other residual systematic effects. Such tests are unique to our 
approach, as the interpretation of analogous active shear estima- 
tors would be hindered by the need to perform a separate PSF de- 
convolution for each estimator, and possibly removing any shared 
defects. 

We have demonstrated the performance of one of our key 



shear estimators via blind tests that were part of the Shear TEst- 
ing Programme (STEP; Massey et al., 2007a). We shall compare 
the practical performance of our remaining shear estimators in a 
future round of STEP simulations. We are also implementing an 
option to input a known flexion signal in the image simulation suite 
of Massey et al. (2004). We are planning a smaller-scale, flexion 
version of STEP, to calibrate the performance of emerging weak 
flexion estimators, including the ones presented in this work as 
well as others presented elsewhere. In the mean time, a complete 
IDL software package capable of performing the shapelet image 
decomposition, including the weak lensing manipulation and anal- 
ysis described in this paper, is available from the shapelets web site 
http : / /www . astro. caltech.edu/~rjm/shapelets. 
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APPENDIX A: REDUCED SHEAR 

Al Idealized case - isophotal or no weighting 

For the purposes of the following discussion we here reproduce 
much of the work of Schneider & Seitz (1995). We use x s and else- 
where the suffix s to denote coordinates and quantities in the galaxy 
source plane, and x and no suffix to denote coordinates and quan- 
tities in the lensed image plane. Let I 3 (x s ) be the surface bright- 
ness distribution of the source and let W(I 3 ) be some weighting 
function of the surface brightness. For the case of no weighting, 
W(I S )=I S . 

We define the centre of the source by 



_ _ J x s W(I s (x s )) dx 3 dy s 
Xs = JW(I B (x s ))dx s dy s ' 

and the quadropole matrix of the source by 

( .) = J(Ax s UAx a )jW(I s (x s )) dx a dy s 
J W(I s {xs)) dx s dy a 



Q 



(Al) 



(A2) 



where Ax s = x s — x s . To describe the shape (including orienta- 
tion) of a source, we define the complex ellipticity, 

_ {Qtf-Q$)+2iQl 



(A3) 



The gravitational imaging of a general source is described by 
the lens equation 



= x — ol(x) 



(A4) 



Where the mass distribution is smooth on scales of galaxy images, 
the imaging of the source can be approximated by the locally lin- 
earized lens mapping 



dxdy = A(x) dx s dy s , 
where 

1 — K — 71 



A(x) = 



-72 



—72 

1 - K + 71 



(A5) 



(A6) 



the Jacobian of the lens equation (A4). 

Now, we can also define analagous second moments Qij for 
the lensed image of a source 

/ AxiAxjW(I(x)) dxdy , 



Qij 



jW(I{x)) dxdy 



(A7) 



Using (A5), and the fact that image surface brightness is conserved 
under gravitational light deflection so that 



W(I s (x s )) = W(I(x)) . 



(A8) 



Using the linearized mapping we may write Axi = Aij(Ax s )j, 
giving the result 



Q\f = A ik AjiQk 



(A9) 



i.e. that Qij transforms as a tensor for a locally linearized map- 
ping. The applicability of this desirable result rests heavily on the 
condition (A8). We may write the Jacobian as 



A(x) = (1-k) 



9i 



-92 



-92 l + 9i 



(A10) 



where we have defined the reduced shear g = 7/(1 — k). The 
transformation between e and e s can then be written as 



e-2g + g 2 e* 



(All) 



1 + | 5 2 |-2Re { 5 e*} ■ 

We see immediately that the transformation between the source and 
image ellipticities e s and e depends solely on the combination g. 

Incidentally, we can continue this calculation one more step 
and obtain, to first order in g, 

e = E s + 2g-2e(Eigi +£252) , (A12) 
which, averaging over a population ensemble, is 

(e) = (ei>+2(l-(e?» ff i-2(eie 2 >92 (A13) 
+ i((e 2 > + 2(1 - {e 2 2 )) 9 2 + 2(eie 2 >ffi) 
= (2 - {e 2 ))g . (A14) 

A2 More general case - weighting by a function of position 

We noted above that the tensor transformation of Qij relies on 
the invariance under lensing transformation of the weighted sur- 
face brightness distribution, a condition that is only satisfied for 
an isophotal weighting function W = W(I). This schema car- 
ries practical difficulties for noisy images, and in general we wish 
to weight objects by multiplying their image by a fixed function 
W(x), such that 



Qij — 



/ Ax l AxjW(x)I(x) dxdy 
J W{x)I{x) dxdy 



(A15) 



It is from weighted moments such as these that weak lensing shear 
is generally measured, and such moments (with gaussian W func- 
tions) are directly equivalent to combinations of the f n 2 and /„n 
polar shapelet coefficients. 

However, we see instantly that the combination I(x)W(x) is 
no longer invariant under the lensing transformation: 



I{x)W(x) + I s (x 3 )W 3 (x s ) 



(A16) 



in general. This prevents us from writing the transformation be- 
tween Qij and Q\ s ^ in the simple form of (A9). The quadrupole 
moments no longer transform as tensors and we must instead write 



Q 



BijkiQki 



(A17) 



where each element of Bijki depends upon 71, 72, tt, and the func- 
tional forms of I and W. Importantly, because each element of 
Bijki will vary with each of these quantities, we cannot therefore 
assume that the transformation between e s and e will depend solely 
on the combination g for an arbitrary weighting function. 

The differences between (A9) and (A17) are generally as- 
sumed to be small for practical weak lensing purposes. Shapelet 
space is convenient for the calculation of elements of Bijki - For a 
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given weighting function, the transformation may be written ap- 
proximately as a power series in 71, 72, k, and the image mo- 
ments/shapelet coefficients f nm . In this way, shapelets provides 
one method for estimating the generalized reduced shear for each 
galaxy image, a complicated function of 71, 72, K, fnm and W in 
each case. 



APPENDIX B: KAISER FLOW 
Bl Population Response 

In §3.4, we obtained expressions for unweighted, unnormalised 
second moments for each galaxy. We constructed an unnormalised 
size measure qo = FR and two unnormalised polarisation com- 
ponents qi = Fei and §2 = FE2, all of which have strong flux 
dependencies. 

We must now find an estimator for the local shear on an image 
given these polarisations. If we were interested in a shear estimate 
for a single galaxy, we might argue that since the lensed quantities 
we measure are related to unlensed quantities by 



q' = q + 2q l3 j l3 
q' a = q a + 2q j a , 

with a,/3= 1, 2, we could use an estimator for the shear 

(la) 



7. 



2(go) 



(Bl) 
(B2) 

(B3) 



However, when we come to combining shear estimators from 
galaxies with different flux and size properties, this approach is not 
adequate. Firstly, it does not give a prescription for how to opti- 
mally combine the estimates from galaxies with very different flux 
and shape properties. Furthermore, it ignores the fact that, under 
shear, some galaxies will flow out of a cell containing galaxies at 
a given flux, qo and ellipticity, while other galaxies will flow in - 
and these flows may not cancel; a weighting scheme should take 
account of this. 

In seeking to address these issues, we closely follow the ap- 
proach offered by Kaiser (2000; section 3.2), although the fact that 
we are dealing with unweighted moments simplifies our analysis. 

We wish to obtain an estimator for the shear that takes into 
account the shear-induced flow of galaxies in the parameter space 
(F,qo,q 2 ), where F is the flux and q 2 = q a q a is an invariant 
measure of the ellipticity amplitude of an object. We will find it 
convenient to describe q a = qq a , with the unit polarisation vec- 
tor given by q a = {cos <f>, sin </>}, i.e. <j> gives the position angle 
of the galaxy. In this case, we can describe a volume element for 
polarisation by qdqd<j), or ^d(q 2 )d4>. 

Let us consider the distribution of galaxies in this parameter 
space. We will represent sheared distributions as primed quantities. 
If the number density in this parameter space is n, then we can 
describe the conservation of galaxies under shear by 



n'(F',q ,q' 2 ,cj>')dF'dq d(q' 2 )dct>' 



(B4) 



n(F,q ,q , <f>)dFdq d(q )d<j> . 

Note that this differs from Kaiser's (2000) analysis in not requiring 
integration over distinct polarisabilities, as these polarisabilities are 
themselves given by qo and q a in our case, due to using unweighted 
moments. 

We now multiply both sides of this equation 
by W(F',q'o,q' 2 )q' a = W(F,qo + 6q ,q 2 + 8q 2 )(q a + Sq a ), 
where W is an arbitrary function, and integrate over all variables. 



This will ultimately allow us to obtain a relation between the aver- 
age polarisation, the distribution of galaxies, and the shear. Initially 
we find 

J dF'dq'od{q' 2 Wn'W{F', q' 0: q' 2 )q' a = (B5) 
J dFdqod(q 2 )d<t>nW(F + SF, qo + Sqo,q 2 + Sq 2 )(q a + Sq a ) . 

This can be simplified by noting that, because of the isotropy of the 
unsheared population, 



dFdqod{q 2 )d4>nW(F,q ,q 2 )q a = . 



(B6) 



/ 

Making a Taylor expansion of the left hand side of equation (B5), 
we obtain 

f dF'dq'od{q' 2 Wn'W{F', q' 0: q' 2 )q' a 

1 



(B7) 



2 / dW dW 2 

dFdq d(q )d(f>n yW8q a + -^—Sq q a + ^ 2 ^ 8(q )q a 

If we note from equations (Bl) and (B2) that Sqo = 2qg^g, 
Sq a = 2go7a, and 8(q 2 ) = ^qojgqg, we can integrate the above 
expression by parts to obtain 



/ 



dF'dq'od(q' 2 )d0'n'W(F', q' , q' 2 )q' a 



(B8) 



7/3 



J dFdqod{q 2 )d(j>w(2nqo8 at 
dn dn 



Since W(F', q' , q' 2 ) = W(F, qo, q 2 ) to first order in the shear, we 
can omit it on both sides (we are free to do this as it is arbitrary), to 
obtain a relation between the mean of the sheared galaxies' polari- 
sations, and the galaxy distribution function, sizes, and shapes 



/ 



nq a dFdq d(q )d<j> = 7/3 J |2rig 5 a(3 (B9) 
_2 9/39« _ 4qimqa 1^2) \ dFd 1o<i(q 2 )d(j> . 



We can usefully write this in terms of an average only over 
position angles of galaxies. If we move to writing expressions in 
terms of the density 



l(F, qo, q 2 ) = j d(f>n(F, q ,q 2 ,(f>) , 



(B10) 



and note that the average over position angles {qpq a } = \q 2 8 al 3, 
then we can write the average of q a over position angle only (i.e. at 
fixed F, qo,q 2 ) as 



{qa)F,q ,q 2 = 7a 



1 dn 2 dn 2 

2 * J "5 ^T~2T*"7 

n oqo n o{q*) 



(BID 



where n is n(F, qo,q 2 ) rather than n(F, qo, q 2 , <j)). It will be 
useful to introduce the susceptibility P, where (q a )F,q , q 2 = 
P{F,qo,q 2 )~l a with 



U/TP 2, „ 1 dn 2 dn 2 

P(F, qo, q ) = 2q ■^T~2\1<- > 1 

n dqo n d(q 2 ) 



(B12) 



We have therefore obtained the appropriate polarisation to use 
as a function of flux, size and shape for an ensemble of galaxies. 
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Hence we can construct a shear estimator for galaxies in a particular 
small cell in (F, qo,q 2 ) space 



-cell 
7a 



— y 

NP ^ 



(B13) 



gals in cell 



where N is the number of galaxies in the cell in question. How- 
ever, we would like an estimator which did not require the split- 
ting of galaxies into cells in parameter space, and which optimally 
combines the estimators from different galaxies. We discuss this 
problem in the next section. 



B2 Optimal Weighting 

Here we will discuss the optimal weighting of shear estimators for 
a spatial cell-average shear. Again, we are following the work of 
Kaiser (2000). 

For our parameter-space cell shear estimate given in equation 
(B13) above, we can find the estimator variance 



N 2p2 



or 



N 2p2 



where the final equality assumes that galaxy polarisations are es- 
sentially uncorrelated in the weak shear regime. Thus we obtain 



((7 CoU ) 2 ) = 



NP 2 



(B15) 



Since parameter-space cells provide shear estimates which are 
uncorrelated from cell to cell, the optimal weighting W C eii is pro- 
portional to l/((7 cc11 ) 2 ) = NP 2 /q 2 , as then the overall estimator 
variance will be minimised. So the final total shear estimate for a 
small spatial aread will be given by 



* total 
7a 



Ecells^WXE 



gals in cell 



q a /NP) 



^-^galaxics Q^ a 



E 

where Q = P/q. 



galaxies 



Q 2 
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APPENDIX C: THE MOST NEARLY LINEAR SHAPE 
ESTIMATOR IN REAL SPACE 

CI Simple polarisation estimators 

The simplest polarisation estimator can be constructed for a galaxy 
image I(x,y) as 



x 2 - V 2 ) I{ X ,V) dxdy 



P2 = j; jj 2xy I(x,y) dxdy 



(CI) 
(C2) 



where the flux F = J J I(x,y) dxdy. The diagonal components 
of the shear susceptibility tensor for this estimator take the simple 
form 2R 2 (evaluated without the weight), and we recover the shear 
estimator ^unweighted from §3.5. The F factor could also have been 
incorporated directly into the shear susceptibility factor; then both 
terms are formally linear, and it is only at the point of forming a 
shear estimator that any ratios need to be taken. 

If we ignore the correction for PSF convolution, the KSB 
method can also be recast in these simple terms. This requires 



Gaussian-weighted quadrupole ellipticities 

Iff - S - ± T~ 

Pi = -j^ // O 2 - y 2 ) e 2r <> I{x,y) 

2 , 2 

Iff ~ -A 

-jp II 2x V e 2r " I(x,y)dxdy, 



dxdy 



P2 = 
where 



R 2 = 



(x 2 + y 2 )e- (x2+y2)/(2r ^I(x,y) dxdy 



(C3) 
(C4) 

(C5) 



and r g is the scale size of a Gaussian weight function. This is in- 
troduced to make sure the integrals converge in a noisy image, and 
to eliminate the effects of neighbouring objects. Unfortunately, this 
weight function complicates the correction for the PSF, and makes 
the corresponding P sh tensor messy (see equations (5-2) to (5-4) 
in KSB). Introducing a ratio of moments at this early stage reduces 
the dynamic range of the ellipticities and the shear susceptibility, 
but also exacerbates the noise. KSB also derive a correction for the 
effects of PSF convolution on a galaxy's shape, with this Gaussian 
and assumptions about the form of the PSF built-in. 

C2 General linear shape estimators 

We can generalise this procedure by defining two arbitrary weight 
functions Wi(x, y), with i € {1, 2}, that can be used to define two 
linear polarisations 



Pi = jj Wi(x,y) I(x,y) dxdy 

The coordinate system is then distorted by a shear 

xl \ = ( 1 + 7i 72 \ / x \ 

y ) V 72 1 - 71 ) \ V J 

and our shape estimators for the observed image become 



Pi = jj Wi(x,y) 



l(x,y) 



(C6) 



(C7) 



(C8) 



dl dl dl dl 

-71^-5- + 7iy « 722;^ 722/ 

ox oy ay 



dx 



dxdy . 



Integrating each term by parts, and including a boundary condition 
on the rapid convergence of the image to zero at large radii, we 
obtain 



- int . 

Pi = Pi + 71 



+ 72 



// 
// 



dm 

dx 



dWi , 

x ~n — h y 

oy 



dW l 

y ^y~ 

dWj 
dx 



I(x,y)dxdy (C9) 
I{x,y) dxdy . 



This pair of integrals, for each of the two weight functions, 
make up the four coefficients in the shear susceptibility tensor. This 
procedure can also be followed in polar coordinates, where we find 



Pi = P. 

+ 



(CIO) 



+ 72 



II 



r sin 26— h cos 26^—- 

Or 30 



I(r,6) rdrdO 
I(r, 6) rdrdO . 



C3 Shear-invariant shear susceptibility 

As was the case in section CI, it is always impossible to form a 
completely linear shear estimator, since bright galaxies would then 



Weak gravitational shear and flexion with polar shapelets 



yield larger shear estimators than faint ones; it will always be nec- 
essary to normalise a shear estimator by something proportional 
to the object's flux. However, we can construct one ellipticity for 
which the shear susceptibility tensor is diagonal and whose diago- 
nal coefficients are exactly equal to that flux (this is equivalent to in- 
cluding a factor of 1/F in an ellipticity estimator that has P 1 = 1). 
This will solve all three problems raised at the beginning of §3, be- 
cause the flux is the most easily measured quantity of an object, 
the matrix inversion can be replaced by a division, and the flux is 
not changed under a shear (nor will the overall shear estimator be 
changed by lensing magnification). 

To achieve this shear estimator with the desired integrals from 
equation (C9), we require 

= J (Cll) 
(C12) 
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dx 
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dW 2 
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dW 2 

x a = 
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dW 2 
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dW! 
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dW x 
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dx 


x 2 + y 2 


dy 


x 2 + y 2 


dW 2 
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dW 2 
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dx 


x 2 + y 2 
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x 2 + y 2 
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or in polar coordinates 

dWi cos (26») dWi 



-sin (26>) (C15) 
cos (26») . (CI 6) 



dr r d8 

dW 2 sin (26 1 ) dW 2 
dr r ' d6 

These equations are inconsistent. Therefore no continuous, analytic 
function exists with all the properties desired for a fully linear shear 
estimator. However, a series approximation that tends to these re- 
quirements is given by the expansion in §3.5. 

This paper has been produced using the Royal Astronomical Soci- 
ety/Blackwell Science MeX style file. 



